Studies 1-4: Figures

Luhrmann, Weisman, et al.

In this notebook, we create figures for Studies 1-4.

source("./scripts_general/dependencies.R")
source("./scripts_general/custom_funs.R")
source("./scripts_general/var_recode_contrast.R")
source("./study1/scripts_s1/s1_var_groups.R")
source("./study2/scripts_s2/s2_var_groups.R")
source("./study3/scripts_s3/s3_var_groups.R")
source("./study4/scripts_s4/s4_var_groups.R")

setwd("./study1/analysis/")

The working directory was changed to /Users/anthrouser/Documents/Projects/Mind and Spirit overview paper/mindspiritquant/study1/analysis inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

source("../../scripts_general/data_load.R")

Parsed with column specification:
cols(
  .default = col_double(),
  study = col_character(),
  date = col_date(format = ""),
  researcher = col_character(),
  country = col_character(),
  site = col_character(),
  religion = col_character(),
  subject_gender = col_character(),
  subject_job = col_character(),
  subject_schedule = col_character(),
  subject_livedhere = col_character(),
  subject_lang = col_character(),
  subject_marital = col_character(),
  subject_hs = col_character(),
  subject_liveswith = col_character(),
  servicesperweek = col_character()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  .default = col_double(),
  study = col_character(),
  date = col_date(format = ""),
  researcher = col_character(),
  country = col_character(),
  site = col_character(),
  religion = col_character(),
  subject_gender = col_character(),
  subject_job = col_character(),
  subject_schedule = col_character(),
  subject_livedhere = col_character(),
  subject_lang = col_character(),
  subject_marital = col_character(),
  subject_hs = col_character(),
  subject_liveswith = col_character(),
  servicesperweek = col_character()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  .default = col_double(),
  ctry = col_character(),
  wher = col_character(),
  recr = col_character(),
  whoc = col_character(),
  demo_chur = col_character(),
  demo_ethn = col_character(),
  demo_maj = col_character(),
  demo_pocc = col_character(),
  demo_rlgn = col_character(),
  demo_sex = col_character()
)
See spec(...) for full column specifications.
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  .default = col_character(),
  X1 = col_double(),
  epi_subj = col_double(),
  epi_demo_age = col_double(),
  epi_demo_ses_num = col_double(),
  epi_demo_howr_num = col_double(),
  epi_demo_affr_num = col_double(),
  epi_demo_tung_num = col_double(),
  epi_demo_ytng_num = col_double(),
  epi_demo_affr_cat = col_logical(),
  epi_demo_tung_cat = col_logical(),
  epi_demo_ytng_cat = col_logical(),
  epi_version = col_double(),
  epi_charc = col_double()
)
See spec(...) for full column specifications.
2 parsing failures.
 row       col expected  actual                           file
1025 epi_charc a double Unclear '../../study3/data/d_demo.csv'
1027 epi_charc a double Unclear '../../study3/data/d_demo.csv'
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  score = col_double()
)
Joining, by = c("epi_ctry", "epi_subj")
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  score = col_double()
)
Joining, by = c("epi_ctry", "epi_subj")
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  question = col_character(),
  response = col_double(),
  order = col_double(),
  question_text = col_character()
)
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  question = col_character(),
  response = col_double(),
  order = col_double(),
  question_text = col_character()
)
Joining, by = c("epi_ctry", "epi_subj", "question", "response", "order", "question_text")
Joining, by = c("epi_ctry", "epi_subj")
Column `epi_ctry` joining character vector and factor, coercing into character vectorJoining, by = "epi_subj"
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  .default = col_double(),
  p7_ctry = col_character(),
  p7_abs_check = col_character(),
  p7_dse_check = col_character(),
  p7_se_check = col_character(),
  p7_unev_check = col_character(),
  p7_exsen_check = col_character(),
  p7_por_check = col_character(),
  p7_mm_check = col_character(),
  p7_dem_sex = col_character(),
  p7_dem_pocc = col_character(),
  p7_dem_major = col_character(),
  p7_dem_ethnicity = col_character(),
  p7_dem_rur.urb = col_character(),
  p7_dem_affrd.basics = col_character(),
  p7_dem_religion = col_character(),
  p7_dem_church = col_character(),
  p7_dem_holy.tung.gif = col_character(),
  p7_abs_child.exp_cat = col_logical(),
  p7_abs_poetic_cat = col_logical(),
  p7_abs_tv.real_cat = col_logical()
  # ... with 162 more columns
)
See spec(...) for full column specifications.
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("study", "p7_ctry", "p7_subj", "abs_score", "cog_score", "ctl_score", "dse_score", "hall_score", "para_score", "por_score", "pv_score", "spev_score")
Joining, by = c("study", "p7_ctry", "p7_subj", "abs_score", "cog_score", "ctl_score", "dse_score", "hall_score", "para_score", "por_score", "pv_score", "spev_score")

d_all <- d1 %>%
  select(study, country, site, religion, subject_id, 
         pv_score, abs_score, spev_score) %>%
  mutate(religion = recode_factor(religion,
                                  "local" = "Local Religion",
                                  "charismatic" = "Charismatic Christianity"),
         # rescale to 0-1
         pv_score = pv_score/3) %>%
  full_join(d2 %>% 
              select(study, country, subj, 
                     abs_score, spev_score, dse_score) %>% 
              rename(subject_id = subj) %>%
              # rescale to 0-1
              mutate(spev_score = spev_score/4,
                     dse_score = dse_score/5,
                     religion = "General Population")) %>%
  full_join(d3 %>% 
              select(study, epi_ctry, epi_sample, epi_subj, 
                     por_score, spirit_score) %>%
              rename(country = epi_ctry, 
                     religion = epi_sample,
                     subject_id = epi_subj,
                     spev_score = spirit_score) %>%
              mutate(religion = recode_factor(religion,
                                              "general population" = "General Population",
                                              "charismatic" = "Charismatic Christianity"))) %>%
  full_join(d4 %>%
              select(study, p7_ctry, p7_subj, 
                     por_score, pv_score, abs_score, spev_score, dse_score) %>%
              rename(country = p7_ctry,
                     subject_id = p7_subj) %>%
              # rescale to 0-1
              mutate(por_score = por_score/2,
                     pv_score = pv_score/3,
                     spev_score = spev_score/4,
                     dse_score = dse_score/5,
                     religion = "General Population")) %>%
  mutate(religion = factor(religion, 
                           levels = c("General Population", 
                                      "Local Religion",
                                      "Charismatic Christianity")),
         study = gsub("study", "Study", study)) 

Joining, by = c("study", "country", "religion", "subject_id", "abs_score", "spev_score")
Column `religion` joining factor and character vector, coercing into character vectorJoining, by = c("study", "country", "religion", "subject_id", "spev_score")
Column `religion` joining character vector and factor, coercing into character vectorJoining, by = c("study", "country", "religion", "subject_id", "pv_score", "abs_score", "spev_score", "dse_score", "por_score")

Summary statistics

By study

Study 4

d_sum_s1 <- d_all %>%
  filter(study == "Study 1") %>%
  group_by(country, religion, site) %>%
  summarise_at(vars(spev_score, pv_score, abs_score),
               funs(mean = mean(., na.rm = T), sd = sd(., na.rm = T))) %>%
  ungroup()

d_rsq_s1 <- data.frame(
  var = c("spev_score", "pv_score", "abs_score"),
  rsq_fix = c(rsq_fun(lm(spev_score ~ country, d_all %>% filter(study == "Study 1"))),
              rsq_fun(lm(pv_score ~ country, d_all %>% filter(study == "Study 1"))),
              rsq_fun(lm(abs_score ~ country, d_all %>% filter(study == "Study 1")))),
  rsq_ran = c(rsq_fun(lmer(spev_score ~ 1 + (1 | country), 
                         d_all %>% filter(study == "Study 1"))),
              rsq_fun(lmer(pv_score ~ 1 + (1 | country), 
                         d_all %>% filter(study == "Study 1"))),
              rsq_fun(lmer(abs_score ~ 1 + (1 | country), 
                         d_all %>% filter(study == "Study 1")))))
d_rsq_s1

plot_s1_spev <- d_all %>% # XX BOOKMARK
  filter(study == "Study 1") %>%
  ggplot(aes(x = country, y = spev_score, 
             color = country, fill = country, 
             shape = site)) +
  facet_grid(~ religion) +
  geom_point(position = position_jitterdodge(jitter.width = 0.25,
                                             jitter.height = 0,
                                             dodge.width = 0.5), 
             alpha = 0.25, show.legend = F) +
  geom_pointrange(data = . %>%
                    group_by(country, religion, site) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  position = position_dodge(width = 0.5),
                  color = "black") +
  geom_text(data = d_sum_s1 %>%
              mutate_at(vars(-country, -religion, -site), 
                        funs(format(round(., 2), nsmall = 2))),
            aes(y = 1, label = paste0(spev_score_mean, "\n(", spev_score_sd, ")")),
                  position = position_dodge(width = 0.75),
            color = "black", size = 2.5, vjust = 1) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_shape_manual(values = 21:24) +
  scale_y_continuous(limits = c(NA, 1), breaks = seq(0, 1, 0.25)) +
  theme(legend.position = "bottom") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  labs(x = "Country", y = "Spiritual Events", shape = "Site")

fig_0_title <- ggdraw() + 
  draw_label("STUDY 1", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 35))

fig_0 <- plot_grid(
  fig_0_title,
  plot_grid(plot_s1_spev, 
            ncol = 1),
  ncol = 1, rel_heights = c(1, 20))

Removed 4 rows containing missing values (geom_point).

# fig_0

fig_0

Study 4

rsq_fun <- function(reg){
  reg_class <- class(reg)
  res <- rsquared(reg)
  
  if (grepl("lmerMod", reg_class)) {
    rsq <- res$Conditional
  } else {
    rsq <- res$R.squared
  }
  return(rsq)
}

d_sum_s4 <- d_all %>%
  filter(study == "Study 4") %>%
  group_by(country) %>%
  summarise_at(vars(spev_score, dse_score, pv_score, por_score, abs_score),
               funs(mean = mean(., na.rm = T), sd = sd(., na.rm = T))) %>%
  ungroup()

d_rsq_s4 <- data.frame(
  var = c("spev_score", "dse_score", "pv_score", "por_score", "abs_score"),
  rsq_fix = c(rsq_fun(lm(spev_score ~ country, d_all %>% filter(study == "Study 4"))),
              rsq_fun(lm(dse_score ~ country, d_all %>% filter(study == "Study 4"))),
              rsq_fun(lm(pv_score ~ country, d_all %>% filter(study == "Study 4"))),
              rsq_fun(lm(por_score ~ country, d_all %>% filter(study == "Study 4"))),
              rsq_fun(lm(abs_score ~ country, d_all %>% filter(study == "Study 4")))),
  rsq_ran = c(rsq_fun(lmer(spev_score ~ 1 + (1 | country), 
                         d_all %>% filter(study == "Study 4"))),
              rsq_fun(lmer(dse_score ~ 1 + (1 | country), 
                         d_all %>% filter(study == "Study 4"))),
              rsq_fun(lmer(pv_score ~ 1 + (1 | country), 
                         d_all %>% filter(study == "Study 4"))),
              rsq_fun(lmer(por_score ~ 1 + (1 | country), 
                         d_all %>% filter(study == "Study 4"))),
              rsq_fun(lmer(abs_score ~ 1 + (1 | country), 
                         d_all %>% filter(study == "Study 4")))))
d_rsq_s4

plot_s4_spev <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = country, y = spev_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  geom_text(data = d_sum_s4 %>%
              mutate_at(vars(-country), funs(format(round(., 2), nsmall = 2))),
            aes(y = 1, label = paste0(spev_score_mean, "\n(", spev_score_sd, ")")),
            color = "black", size = 2.5, vjust = 1) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Spiritual Events")

plot_s4_dse <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = country, y = dse_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(dse_score, na.rm = T),
                              sd = sd(dse_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  geom_text(data = d_sum_s4 %>%
              mutate_at(vars(-country), funs(format(round(., 2), nsmall = 2))),
            aes(y = 1, label = paste0(dse_score_mean, "\n(", dse_score_sd, ")")),
            color = "black", size = 2.5, vjust = 1) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Daily Spiritual Experiences")

plot_s4_pv <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = country, y = pv_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(pv_score, na.rm = T),
                              sd = sd(pv_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  geom_text(data = d_sum_s4 %>%
              mutate_at(vars(-country), funs(format(round(., 2), nsmall = 2))),
            aes(y = 1, label = paste0(pv_score_mean, "\n(", pv_score_sd, ")")),
            color = "black", size = 2.5, vjust = 1) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Porosity Vignettes")

plot_s4_por <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = country, y = por_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(por_score, na.rm = T),
                              sd = sd(por_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  geom_text(data = d_sum_s4 %>%
              mutate_at(vars(-country), funs(format(round(., 2), nsmall = 2))),
            aes(y = 1, label = paste0(por_score_mean, "\n(", por_score_sd, ")")),
            color = "black", size = 2.5, vjust = 1) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Porosity Scale")

plot_s4_abs <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = country, y = abs_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(abs_score, na.rm = T),
                              sd = sd(abs_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  geom_text(data = d_sum_s4 %>%
              mutate_at(vars(-country), funs(format(round(., 2), nsmall = 2))),
            aes(y = 1, label = paste0(abs_score_mean, "\n(", abs_score_sd, ")")),
            color = "black", size = 2.5, vjust = 1) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Absorption")

fig_2_title <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 35))

fig_2 <- plot_grid(
  fig_2_title,
  plot_grid(plot_s4_spev, plot_s4_dse, 
            plot_s4_pv, plot_s4_por, 
            plot_s4_abs, NULL, 
            ncol = 2, labels = c("A", "B", "C", "D", "E")),
  ncol = 1, rel_heights = c(1, 20))

Removed 2 rows containing missing values (geom_point).

# fig_2

fig_2

Spiritual Events scores

Study 1

d1 %>%
  ggplot(aes(x = country, y = spev_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Study 1: Spiritual Events score (range: 0-1)",
       caption = "Error bars are ±1 standard deviation from the mean")

d_all %>%
  filter(study == "Study 1") %>%
  ggplot(aes(x = country, y = spev_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country, religion) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd,
                      shape = religion),
                  color = "black",
                  show.legend = T) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_shape_manual(values = 21:24) +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.25)) +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  labs(x = "Country", y = "Study 1: Spiritual Events score (range: 0-1)",
       caption = "Error bars are ±1 standard deviation from the mean")

Study 2

d2 %>%
  ggplot(aes(x = country, y = spev_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 4), breaks = seq(0, 4, 1)) +
  labs(x = "Country", y = "Study 2: Spiritual Events score (range: 0-4)",
       caption = "Error bars are ±1 standard deviation from the mean")

Study 3

d3 %>%
  # correct for scaling in original dataset
  mutate(spev_score = spirit_score * 4) %>%
  ggplot(aes(x = epi_ctry, y = spev_score, color = epi_ctry, fill = epi_ctry)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(epi_ctry) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 4), breaks = seq(0, 4, 1)) +
  labs(x = "Country", y = "Study 3: Spiritual Events score (range: 0-4)",
       caption = "Error bars are ±1 standard deviation from the mean")

Study 4

d4 %>%
  ggplot(aes(x = p7_ctry, y = spev_score, color = p7_ctry, fill = p7_ctry)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(p7_ctry) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 4), breaks = seq(0, 4, 1)) +
  labs(x = "Country", y = "Study 4: Spiritual Events score (range: 0-4)",
       caption = "Error bars are ±1 standard deviation from the mean")

All studies

d_all %>%
  ggplot(aes(x = country, y = spev_score, color = country, fill = country,
             group = study)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.25,
                                             jitter.height = 0,
                                             dodge.width = 0.75), 
             alpha = 0.1) +
  geom_pointrange(data = d_all %>%
                    group_by(study, country) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, shape = study), 
                  position = position_dodge(width = 0.75),
                  color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom") +
  labs(x = "Country", y = "Spiritual Events score (rescaled to 0-1)",
       shape = "Study", 
       caption = "Error bars are ±1 standard deviation from the mean")

d_spev_sum <- d_all %>%
  filter(!is.na(spev_score)) %>%
  group_by(study, country, religion) %>%
  summarise(mean = mean(spev_score, na.rm = T),
            sd = sd(spev_score, na.rm = T),
            n = n()) %>%
  ungroup()

d_all %>%
  ggplot(aes(x = country, y = spev_score, 
             color = country, fill = country,
             group = religion)) +
  facet_grid(~ study, scales = "free", space = "free") +
  geom_point(position = position_jitterdodge(jitter.width = 0.8,
                                             jitter.height = 0.02,
                                             dodge.width = 0.75), 
             alpha = 0.15) +
  geom_pointrange(data = d_spev_sum,
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, 
                      shape = religion), 
                  position = position_dodge(width = 0.75),
                  fill = "black",
                  color = "black") +
  geom_text(data = d_spev_sum %>%
              mutate(ypos = case_when(
                grepl("charismatic", tolower(religion)) ~ mean + sd + 0.05,
                TRUE ~ mean - sd - 0.05)),
            aes(y = ypos, label = paste0("n=", n)), 
            position = position_dodge(width = 0.75),
            size = 3, color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "Country", y = "Spiritual Events score (rescaled to 0-1)",
       # caption = "Error bars are ±1 standard deviation from the mean",
       shape = "Religion")

d_all %>%
  filter(study == "Study 1") %>%
  ggplot(aes(x = country, y = spev_score, 
             color = country, fill = country,
             group = religion)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.8,
                                             jitter.height = 0.02,
                                             dodge.width = 0.75), 
             alpha = 0.15) +
  geom_pointrange(data = d_spev_sum %>%
                    filter(study == "Study 1"),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, 
                      shape = religion), 
                  position = position_dodge(width = 0.75),
                  fill = "black",
                  color = "black") +
  geom_text(data = d_spev_sum %>%
              filter(study == "Study 1") %>%
              mutate(ypos = case_when(
                grepl("charismatic", tolower(religion)) ~ mean + sd + 0.05,
                TRUE ~ mean - sd - 0.05)),
            aes(y = ypos, label = paste0("n=", n)), 
            position = position_dodge(width = 0.75),
            size = 3, color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom") +
  labs(x = "Country", y = "Spiritual Events score (rescaled to 0-1)",
       # caption = "Error bars are ±1 standard deviation from the mean",
       shape = "Religion")

Porosity scores

d_por_sum <- d_all %>%
  gather(por_scale, por_score, c(pv_score, por_score)) %>%
  mutate(por_scale = recode(por_scale,
                            "pv_score" = "Porosity Vignettes",
                            "por_score" = "Porosity Scale")) %>%
  filter(!is.na(por_score)) %>%
  group_by(study, country, religion, por_scale) %>%
  summarise(mean = mean(por_score, na.rm = T),
            sd = sd(por_score, na.rm = T),
            n = n()) %>%
  ungroup()

Factor `religion` contains implicit NA, consider using `forcats::fct_explicit_na`

d_all %>% 
  gather(por_scale, por_score, c(pv_score, por_score)) %>%
  mutate(por_scale = recode(por_scale,
                            "pv_score" = "Porosity Vignettes",
                            "por_score" = "Porosity Scale")) %>%
  filter(!is.na(por_score)) %>%
  ggplot(aes(x = country, y = por_score, 
             color = country, fill = country,
             group = religion)) +
  facet_grid(cols = vars(study, por_scale), scales = "free", space = "free") +
  geom_point(position = position_jitterdodge(jitter.width = 0.8,
                                             jitter.height = 0.02,
                                             dodge.width = 0.75), 
             alpha = 0.15) +
  geom_pointrange(data = d_por_sum,
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, 
                      shape = religion), 
                  position = position_dodge(width = 0.75),
                  fill = "black",
                  color = "black") +
  geom_text(data = d_por_sum %>%
              mutate(ypos = case_when(
                grepl("charismatic", tolower(religion)) ~ mean + sd + 0.05,
                TRUE ~ mean - sd - 0.05)),
            aes(y = ypos, label = paste0("n=", n)), 
            position = position_dodge(width = 0.75),
            size = 3, color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "Country", y = "Porosity score (rescaled to 0-1)",
       # caption = "Error bars are ±1 standard deviation from the mean",
       shape = "Religion")

Absorption scores

d_abs_sum <- d_all %>%
  filter(!is.na(abs_score)) %>%
  group_by(study, country, religion) %>%
  summarise(mean = mean(abs_score, na.rm = T),
            sd = sd(abs_score, na.rm = T),
            n = n()) %>%
  ungroup()

Factor `religion` contains implicit NA, consider using `forcats::fct_explicit_na`

d_all %>% 
  filter(!is.na(abs_score)) %>%
  ggplot(aes(x = country, y = abs_score, 
             color = country, fill = country,
             group = religion)) +
  facet_grid(. ~ study, scales = "free", space = "free") +
  geom_point(position = position_jitterdodge(jitter.width = 0.8,
                                             jitter.height = 0.02,
                                             dodge.width = 0.75), 
             alpha = 0.15) +
  geom_pointrange(data = d_abs_sum,
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, 
                      shape = religion), 
                  position = position_dodge(width = 0.75),
                  fill = "black",
                  color = "black") +
  geom_text(data = d_abs_sum %>%
              mutate(ypos = case_when(
                grepl("charismatic", tolower(religion)) ~ mean + sd + 0.05,
                TRUE ~ mean - sd - 0.05)),
            aes(y = ypos, label = paste0("n=", n)), 
            position = position_dodge(width = 0.75),
            size = 3, color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "Country", y = "Absorption score (rescaled to 0-1)",
       # caption = "Error bars are ±1 standard deviation from the mean",
       shape = "Religion")

r1_spev <- lm(spev_score_std ~ country, d1)
r1_pv <- lm(pv_score_std ~ country, d1)
r1_abs <- lm(abs_score_std ~ country, d1)

r2_spev <- lm(spev_score_std ~ country, d2)
r2_dse <- lm(dse_score_std ~ country, d2)
r2_abs <- lm(abs_score_std ~ country, d2)

r3_spev <- lm(spirit_score_std ~ epi_ctry, d3)
r3_por <- lm(por_score_std ~ epi_ctry, d3)

r4_spev <- lm(spev_score_std ~ p7_ctry, d4)
r4_dse <- lm(dse_score_std ~ p7_ctry, d4)
r4_por <- lm(por_score_std ~ p7_ctry, d4)
r4_pv <- lm(pv_score_std ~ p7_ctry, d4)
r4_abs <- lm(abs_score_std ~ p7_ctry, d4)

df1 <- data.frame(study = "study 1",
                  var = c("spiritual experience", "porosity", "absorption"),
                  scale = c("spiritual events", "porosity vignettes", "absorption"),
                  rsq = c(rsquared(r1_spev)$R.squared, 
                          rsquared(r1_pv)$R.squared, 
                          rsquared(r1_abs)$R.squared))

df2 <- data.frame(study = "study 3",
                  var = c("spiritual experience", "spiritual experience", "absorption"),
                  scale = c("spiritual events", "DSE", "absorption"),
                  rsq = c(rsquared(r2_spev)$R.squared, 
                          rsquared(r2_dse)$R.squared, 
                          rsquared(r2_abs)$R.squared))

df3 <- data.frame(study = "study 2",
                  var = c("spiritual experience", "porosity"),
                  scale = c("spiritual events", "porosity scale"),
                  rsq = c(rsquared(r3_spev)$R.squared, 
                          rsquared(r3_por)$R.squared))

df4 <- data.frame(study = "study 4",
                  var = c("spiritual experience", "spiritual experience", 
                          "porosity", "porosity", "absorption"),
                  scale = c("spiritual events", "DSE", 
                            "porosity scale", "porosity vignettes", "absorption"),
                  rsq = c(rsquared(r4_spev)$R.squared, 
                          rsquared(r4_dse)$R.squared,
                          rsquared(r4_por)$R.squared, 
                          rsquared(r4_pv)$R.squared, 
                          rsquared(r4_abs)$R.squared))

df_all <- full_join(df1, df2) %>% full_join(df3) %>% full_join(df4) %>%
  mutate(var = factor(var, 
                      levels = c("spiritual experience", 
                                 "porosity", "absorption"))) %>%
  select(var, scale, study, rsq) %>%
  arrange(var, scale, study)

Joining, by = c("study", "var", "scale", "rsq")
Column `study` joining factors with different levels, coercing to character vectorColumn `var` joining factors with different levels, coercing to character vectorColumn `scale` joining factors with different levels, coercing to character vectorJoining, by = c("study", "var", "scale", "rsq")
Column `study` joining character vector and factor, coercing into character vectorColumn `var` joining character vector and factor, coercing into character vectorColumn `scale` joining character vector and factor, coercing into character vectorJoining, by = c("study", "var", "scale", "rsq")
Column `study` joining character vector and factor, coercing into character vectorColumn `var` joining character vector and factor, coercing into character vectorColumn `scale` joining character vector and factor, coercing into character vector

df_all %>% 
  mutate(percent_exp = paste0(round(rsq * 100), "%")) %>%
  select(-rsq) %>%
  spread(study, percent_exp) %>%
  mutate_at(vars(starts_with("study")),
            funs(case_when(is.na(.) ~ ".", 
                           TRUE ~ .))) %>%
  kable() %>% 
  kable_styling() %>%
  collapse_rows(1:3)

var	scale	study 1	study 2	study 3	study 4
spiritual experience	DSE	.	.	50%	50%
	spiritual events	22%	24%	38%	31%
porosity	porosity scale	.	61%	.	39%
	porosity vignettes	37%	.	.	27%
absorption	absorption	10%	.	11%	9%

r1_spev_pv <- lm(spev_score_std ~ pv_score_std, d1)
r1_spev_abs <- lm(spev_score_std ~ abs_score_std, d1)

r2_spev_abs <- lm(spev_score_std ~ abs_score_std, d2)
r2_dse_abs <- lm(dse_score_std ~ abs_score_std, d2)

r3_spev_por <- lm(spirit_score_std ~ por_score_std, d3)

r4_spev_por <- lm(spev_score_std ~ por_score_std, d4)
r4_dse_por <- lm(dse_score_std ~ por_score_std, d4)
r4_spev_pv <- lm(spev_score_std ~ pv_score_std, d4)
r4_dse_pv <- lm(dse_score_std ~ pv_score_std, d4)
r4_spev_abs <- lm(spev_score_std ~ abs_score_std, d4)
r4_dse_abs <- lm(dse_score_std ~ abs_score_std, d4)

df1b <- data.frame(study = "study 1",
                   outcome = "spiritual events",
                   predictor = c("porosity vignettes", "absorption"),
                   rsq = c(rsquared(r1_spev_pv)$R.squared, 
                           rsquared(r1_spev_abs)$R.squared))

df2b <- data.frame(study = "study 3",
                   outcome = c("spiritual events", "DSE"),
                   predictor = "absorption",
                   rsq = c(rsquared(r2_spev_abs)$R.squared, 
                           rsquared(r2_dse_abs)$R.squared))

df3b <- data.frame(study = "study 2",
                   outcome = "spiritual events",
                   predictor = "porosity scale",
                   rsq = c(rsquared(r3_spev_por)$R.squared))

df4b <- data.frame(study = "study 4",
                   outcome = rep(c("spiritual events", "DSE"), 3),
                   predictor = c(rep("porosity scale", 2), 
                                 rep("porosity vignettes", 2),
                                 rep("absorption", 2)),
                   rsq = c(rsquared(r4_spev_por)$R.squared, 
                           rsquared(r4_dse_por)$R.squared,
                           rsquared(r4_spev_pv)$R.squared, 
                           rsquared(r4_dse_pv)$R.squared, 
                           rsquared(r4_spev_abs)$R.squared,
                           rsquared(r4_dse_abs)$R.squared))

df_allb <- full_join(df1b, df2b) %>% full_join(df3b) %>% full_join(df4b) %>%
  mutate(outcome = factor(outcome,
                          levels = c("spiritual events", "DSE")),
         predictor = factor(predictor,
                            levels = c("porosity vignettes", "porosity scale", "absorption"))) %>%
  select(predictor, outcome, study, rsq) %>%
  arrange(predictor, outcome, study)

Joining, by = c("study", "outcome", "predictor", "rsq")
Column `study` joining factors with different levels, coercing to character vectorColumn `outcome` joining factors with different levels, coercing to character vectorColumn `predictor` joining factors with different levels, coercing to character vectorJoining, by = c("study", "outcome", "predictor", "rsq")
Column `study` joining character vector and factor, coercing into character vectorColumn `outcome` joining character vector and factor, coercing into character vectorColumn `predictor` joining character vector and factor, coercing into character vectorJoining, by = c("study", "outcome", "predictor", "rsq")
Column `study` joining character vector and factor, coercing into character vectorColumn `outcome` joining character vector and factor, coercing into character vectorColumn `predictor` joining character vector and factor, coercing into character vector

df_allb %>% 
  full_join(df_all %>% 
              filter(var == "spiritual experience") %>%
              select(-var) %>%
              rename(outcome = scale) %>%
              mutate(predictor = "country")) %>%
  mutate(predictor = factor(predictor,
                            levels = c("country", "porosity vignettes",
                                       "porosity scale", "absorption"))) %>%
  mutate(percent_exp = paste0(round(rsq * 100), "%")) %>%
  select(-rsq) %>%
  spread(study, percent_exp) %>%
  mutate_at(vars(starts_with("study")),
            funs(case_when(is.na(.) ~ ".", 
                           TRUE ~ .))) %>%
  kable() %>% 
  kable_styling() %>%
  collapse_rows(1:2)

Joining, by = c("predictor", "outcome", "study", "rsq")
Column `predictor` joining factor and character vector, coercing into character vectorColumn `outcome` joining factor and character vector, coercing into character vector

predictor	outcome	study 1	study 2	study 3	study 4
country	DSE	.	.	50%	50%
	spiritual events	22%	24%	38%	31%
porosity vignettes	DSE	.	.	.	28%
	spiritual events	13%	.	.	29%
porosity scale	DSE	.	.	.	44%
	spiritual events	.	39%	.	33%
absorption	DSE	.	.	3%	2%
	spiritual events	11%	.	12%	8%

Relationships

All studies, multipart plot

Spiritual Events only, cowplot

fig_s1_por <- d_all %>%
  filter(study == "Study 1") %>%
  ggplot(aes(x = pv_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Vignettes",
       y = "Spiritual Events",
       color = "Country")
# fig_s1_por

fig_s1_abs <- d_all %>%
  filter(study == "Study 1") %>%
  ggplot(aes(x = abs_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Spiritual Events",
       color = "Country")
# fig_s1_abs

fig_s1_title <- ggdraw() + 
  draw_label("STUDY 1", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 105))

fig_s1 <- plot_grid(
  fig_s1_title,
  plot_grid(fig_s1_por, fig_s1_abs, ncol = 1, labels = c("A", "B")),
  ncol = 1, rel_heights = c(1, 10))

Removed 20 rows containing non-finite values (stat_smooth).Removed 20 rows containing non-finite values (stat_smooth).Removed 20 rows containing missing values (geom_point).Removed 28 rows containing non-finite values (stat_smooth).Removed 28 rows containing non-finite values (stat_smooth).Removed 28 rows containing missing values (geom_point).

# fig_s1

fig_s1_title_vert <- ggdraw() + 
  draw_label("STUDY 1", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s1_vert <- plot_grid(
  fig_s1_title_vert,
  plot_grid(fig_s1_por, fig_s1_abs, ncol = 2, labels = c("A", "B")),
  ncol = 1, rel_heights = c(1, 10))

Removed 20 rows containing non-finite values (stat_smooth).Removed 20 rows containing non-finite values (stat_smooth).Removed 20 rows containing missing values (geom_point).Removed 28 rows containing non-finite values (stat_smooth).Removed 28 rows containing non-finite values (stat_smooth).Removed 28 rows containing missing values (geom_point).

# fig_s1_vert

fig_s2_abs <- d_all %>%
  filter(study == "Study 2") %>%
  ggplot(aes(x = abs_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Spiritual Events",
       color = "Country")
# fig_s2_abs

fig_s2_title <- ggdraw() + 
  draw_label("STUDY 3", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 105))

fig_s2 <- plot_grid(
  fig_s2_title,
  plot_grid(NULL, fig_s2_abs, ncol = 1, labels = c("", "C")),
  ncol = 1, rel_heights = c(1, 10))

Removed 10 rows containing missing values (geom_smooth).

# fig_s2

fig_s2_title_vert <- ggdraw() + 
  draw_label("STUDY 3", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s2_vert <- plot_grid(
  fig_s2_title_vert,
  plot_grid(NULL, fig_s2_abs, ncol = 2, labels = c("", "C")),
  ncol = 1, rel_heights = c(1, 10))

Removed 10 rows containing missing values (geom_smooth).

# fig_s2_vert

fig_s3_por <- d_all %>%
  filter(study == "Study 3") %>%
  ggplot(aes(x = por_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Scale",
       y = "Spiritual Events",
       color = "Country")
# fig_s3_por

fig_s3_title <- ggdraw() + 
  draw_label("STUDY 2", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 105))

fig_s3 <- plot_grid(
  fig_s3_title,
  plot_grid(fig_s3_por, NULL, ncol = 1, labels = c("D", "")),
  ncol = 1, rel_heights = c(1, 10))
# fig_s3

fig_s3_title_vert <- ggdraw() + 
  draw_label("STUDY 2", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s3_vert <- plot_grid(
  fig_s3_title_vert,
  plot_grid(fig_s3_por, NULL, ncol = 2, labels = c("D", "")),
  ncol = 1, rel_heights = c(1, 10))
# fig_s3_vert

fig_s32_vert <- plot_grid(
  plot_grid(
    fig_s3_title_vert, 
    plot_grid(fig_s3_por, labels = c("C")), 
    ncol = 1, rel_heights = c(1, 10)),
  plot_grid(
    fig_s2_title_vert,
    plot_grid(fig_s2_abs, labels = c("D")), 
    ncol = 1, rel_heights = c(1, 10))
)

Removed 10 rows containing missing values (geom_smooth).

fig_s4_por1 <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = pv_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Vignettes",
       y = "Spiritual Events",
       color = "Country")
# fig_s4_por1

fig_s4_por2 <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = por_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Scale",
       y = "Spiritual Events",
       color = "Country")
# fig_s4_por2

fig_s4_abs <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = abs_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Spiritual Events",
       color = "Country")
# fig_s4_abs

fig_s4_title <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 190))

fig_s4 <- plot_grid(
  fig_s4_title,
  plot_grid(plot_grid(fig_s4_por1, fig_s4_por2, ncol = 2, labels = c("E", "F")), 
            plot_grid(NULL, fig_s4_abs, NULL, ncol = 3, rel_widths = c(1, 2, 1), labels = c("", "G", "")), 
            ncol = 1),
  ncol = 1, rel_heights = c(1, 10))

Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 11 rows containing missing values (geom_smooth).

# fig_s4

fig_s4_title_vert <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s4_vert <- plot_grid(
  fig_s4_title_vert,
  plot_grid(plot_grid(fig_s4_por1, fig_s4_por2, ncol = 1, labels = c("E", "F")), 
            plot_grid(NULL, fig_s4_abs, NULL, ncol = 1, rel_heights = c(1, 2, 1), labels = c("", "G", "")), 
            ncol = 2),
  ncol = 1, rel_heights = c(1, 20))

Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 11 rows containing missing values (geom_smooth).

# fig_s4_vert

fig_legend <- get_legend(fig_s4_por1 + theme(legend.position = "bottom"))

Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).

fig_all <- plot_grid(fig_s1, fig_s2, fig_s3, fig_s4, ncol = 4,
                     rel_widths = c(1, 1, 1, 2), scale = 0.95)
fig_all

fig_all_vert <- plot_grid(fig_s1_vert, fig_s32_vert, fig_s4_vert, fig_legend,
                          ncol = 1, rel_heights = c(1, 1, 2, 0.2))
fig_all_vert

Daily Spiritual Experiences only, cowplot

fig_s2_abs <- d_all %>%
  filter(study == "Study 2") %>%
  ggplot(aes(x = abs_score, y = dse_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Daily Spiritual Experiences",
       color = "Country")
# fig_s2_abs

fig_s2_title_vert <- ggdraw() + 
  draw_label("STUDY 3", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s2_vert <- plot_grid(
  fig_s2_title_vert,
  plot_grid(fig_s2_abs, ncol = 1, labels = "B"),
  ncol = 1, rel_heights = c(1, 10))
# fig_s2_vert

fig_s4_por1 <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = pv_score, y = dse_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Vignettes",
       y = "Daily Spiritual Experiences",
       color = "Country")
# fig_s4_por1

fig_s4_por2 <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = por_score, y = dse_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Scale",
       y = "Daily Spiritual Experiences",
       color = "Country")
# fig_s4_por2

fig_s4_abs <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = abs_score, y = dse_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Daily Spiritual Experiences",
       color = "Country")
# fig_s4_abs

fig_s4_title_vert <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s4_por1_vert <- plot_grid(
  fig_s4_title_vert,
  plot_grid(fig_s4_por1, ncol = 1, labels = "A"),
  ncol = 1, rel_heights = c(1, 10))
# fig_s4_por1_vert

fig_s4_por2_vert <- plot_grid(
  fig_s4_title_vert,
  plot_grid(fig_s4_por2, ncol = 1, labels = "C"),
  ncol = 1, rel_heights = c(1, 10))
# fig_s4_por2_vert

fig_s4_abs_vert <- plot_grid(
  fig_s4_title_vert,
  plot_grid(fig_s4_abs, ncol = 1, labels = "D"),
  ncol = 1, rel_heights = c(1, 10))
# fig_s4_abs_vert

fig_legend <- get_legend(fig_s4_por1 + theme(legend.position = "bottom"))

fig_all_vert <- plot_grid(
  plot_grid(fig_s4_por1_vert, fig_s2_vert, 
            fig_s4_por2_vert, fig_s4_abs_vert,
            ncol = 2),
  fig_legend,
  ncol = 1, rel_heights = c(2, 0.2))
fig_all_vert

Other versions

Spiritual Events only, one grid, new layout

d_all %>%
  gather(spirit_scale, spev_score, c(spev_score, dse_score)) %>%
  gather(pred_scale, pred_score, c(por_score, pv_score, abs_score)) %>%
  mutate(spirit_scale = recode_factor(spirit_scale,
                                      "spev_score" = "Spiritual Events",
                                      "dse_score" = "Daily Spiritual Experiences"),
         pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption"),
         pred_type = case_when(pred_scale == "Absorption" ~ "Absorption",
                               grepl("Porosity", pred_scale) ~ "Porosity",
                               TRUE ~ NA_character_),
         pred_type = factor(pred_type, levels = c("Porosity", "Absorption")),
         study_scale = paste(study, pred_scale, sep = ": "),
         study_scale2 = case_when(
           study == "Study 4"  & pred_scale != "Absorption" ~ 
             paste(study, pred_scale, sep = ": "),
           TRUE ~ study),
         study_scale3 = case_when(
           study == "Study 4" & pred_scale == "Porosity Scale" ~ "Porosity Scale",
           study %in% c("Study 1", "Study 4") ~ "Porosity Vignettes",
           study == "Study 2" ~ "Porosity Scale", 
           TRUE ~ " "),
         study_scale3 = factor(study_scale3,
                               levels = c("Porosity Vignettes",
                                          "Porosity Scale", " "))) %>%
  filter(!is.na(pred_score),
         spirit_scale == "Spiritual Events") %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(rows = vars(pred_type), cols = vars(study, study_scale3)) +
  # facet_grid(pred_type ~ study_scale3) +
  geom_point(data = . %>% distinct(study, study_scale3, country, 
                                   pred_type, pred_scale, pred_score,
                                   spirit_scale, spev_score),
             aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Score on predictor scale (Porosity Vignettes, Porosity Scale, or Absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

Spiritual Events only, by study

fig_s1 <- d_all %>%
  filter(study == "Study 1") %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

fig_s2 <- d_all %>%
  filter(study == "Study 2") %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

fig_s3 <- d_all %>%
  filter(study == "Study 3") %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

fig_s4 <- d_all %>%
  filter(study == "Study 4") %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

Spiritual Events only, full grid

d_all %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

Broken down by predictor type and spiritual scale

d_all %>%
  gather(spirit_scale, spev_score, c(spev_score, dse_score)) %>%
  gather(poros_scale, poros_score, c(por_score, pv_score)) %>%
  mutate(spirit_scale = recode_factor(spirit_scale,
                                      "spev_score" = "Spiritual Events",
                                      "dse_score" = "Daily Spiritual Experiences"),
         poros_scale = recode_factor(poros_scale,
                                     "pv_score" = "Porosity Vignettes",
                                     "por_score" = "Porosity Scale"),
         study_scale = paste(study, poros_scale, sep = ": ")) %>%
  filter(!is.na(poros_score)) %>%
  ggplot(aes(x = poros_score, y = spev_score)) +
  facet_grid(spirit_scale ~ study_scale) +
  geom_point(data = . %>% distinct(study, study_scale, country, 
                                   poros_scale, poros_score,
                                   spirit_scale, spev_score),
             aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Score on porosity measure (rescaled to 0-1)",
       y = "Score on spiritual experience measure (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

d_all %>%
  gather(spirit_scale, spev_score, c(spev_score, dse_score)) %>%
  mutate(spirit_scale = recode_factor(spirit_scale,
                                      "spev_score" = "Spiritual Events",
                                      "dse_score" = "Daily Spiritual Experiences"),
         study_scale = paste(study, "Absorption scale", sep = ": ")) %>%
  filter(!is.na(abs_score)) %>%
  ggplot(aes(x = abs_score, y = spev_score)) +
  facet_grid(spirit_scale ~ study_scale) +
  geom_point(data = . %>% distinct(study, study_scale, country, 
                                   abs_score,
                                   spirit_scale, spev_score),
             aes(color = country), alpha = 0.2) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Score on absorption measure (rescaled to 0-1)",
       y = "Score on spiritual experience measure (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")